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Douglas Darling and Barbara Sakowski 
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Nomenclature 

A flow area 

F s momentum source term used in core flow 

equations 

M Mach number 

M R niass source term used in core flow equations 
P local static pressure 

Q s heat source term used in core flow equations 
r c radius of cowl 

r s radius of centerbody 

R t radius of the cowl lip. used to nondimen- 
sionalize the core flow equations 
Re local Reynolds Number per unit length 

R^ radius of the wall (centerbody or cowl), a 
function of x but not a function of y 
RVB coefficient in front of derivatives of advective 
terms in y-direction 

u local velocity component parallel to the core 

flow 

v local velocity component perpendicular to the 

core flow 

T local static temperature 

t time 

x distance, direction of core flow 

y distance from wall, direction perpendicular to 

the core flow 

y 1 distance from wall, transformed by core flow 
stream area 
p density 

non-dimensional variable in core flow equ- 
ations 

non-dimensional variable in boundary layer 
equations 

„ freestream condition, upstream of cowl and 

centerbody 

c core flow, boundary layer edge, condition 

max core A° w value at point with the highest 

Mach number, these values are used to non- 
diniensionalize the boundary layer equations. 

0 stagnation quantity 

Kt reference values used to non-dimensionalize 

core flow equations 


I. Introduction 

Bleeds, bypasses, centerbody translations, and active 
cooling are some of the complexities of modeling flow 
within a supersonic jet engine inlet. When developing 
control laws for inlets a large number of geometries and 
flow conditions need to be modeled. When a large 
number of geometries and flow are to be modeled the 
speed of the modeling technique is very important. It 
would be beneficial to model the bulk of the flow with 
a quick inviscid model. In fact, much of the flow inside 
the inlet of a supersonic jet engine can indeed be treated 
as inviscid. However, at times, the boundary layer does 
have a significant effect on the overall flow within an 
inlet. Thus, we would like a quick and robust method 
to model the boundary layer effects. It was not the 
intent of this work to develop a detailed model of 
shock-wave boundary layer interactions or a detailed 
flow model, which would be required for a detailed inlet 
design. The purpose of this work was to develop a fast 
method to model inlets, including boundary layer ef- 
fects, when the bulk flow properties and overall flow re- 
sponses are of the biggest concern, such as in a controls 
application. 

The Large Perturbation Inlet computer code (LAPIN) 
was developed to model inlets of supersonic jet engines. 
LAPIN modeled the flow within an inlet as one-dimen- 
sional and inviscid. For the most part, this inviscid 
model did a good job of predicting the flow within a 
supersonic inlet. However, in regions of strong adverse 
pressure gradients, near a shock wave, or near the throat 
of the inlet, the inviscid model did not predict the inlet 
flow well. These were the regions where boundary 
layer effects had a significant effect on the overall flow 
within an inlet. 

Roach, et al. [I] showed that an uncoupled, finite differ- 
ence boundary layer scheme is a fast and accurate algo- 
rithm to apply when modeling unsteady, supersonic 
flows. Using this uncoupled scheme we were able to 
model attached or separated and turbulent or laminar 
flows. 

In this paper we will look at the problem of interfac- 
ing this boundary layer model with an inviscid, core 
flow model (LAPIN). Then the results of the combined 



viscous/inviscid model will be discussed and compared 
with results of experiments that had been run on the 
NASA - Lewis 40/60 inlet. 

II. Core Flow and Boundary Layer Models 

The core flow and boundary layer regions of the flow 
field within an inlet were modeled in different manners. 
This was done due to differences in the governing equa- 
tions to model the two regions and because they were 
developed independently. To understand how to inter- 
face the two regions (viscous and inviscid) we will first 
look at each model separately. 

The core flow equations were solved by the program 
LAPIN in conservation form for the flux variables. 

Then the primitive variables were backed out of the flux 
terms. The equations were one-dimensional, inviscid, 
unsteady, and included source terms in each equation. 

The core flow governing equations, used in LAPIN, can 
be written in non-dimensional vector form as— 
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and the nondimensional variables were defined as— 


u r , f t 


(3) 


The reference length. R ( . was the radius of the cowl lip. 
All other reference conditions were taken from the 
upstream sonic conditions. LAPIN solved these equa- 
tions by marching in time, using any of the follow ing 
algorithms: Beam-Warming, hybrid Beam Warming. 
MacCormack’s method, or split characteristics MacCor- 
mack's method. The split characteristics method tended 
to yield the best results and allowed the largest time 
step while maintaining numerical stability. LAPIN was 
set up to model centerbody translations, bleeds, bypass- 
es, heat sources, and mass injection. A detailed descrip- 
tion of the core flow equations and solution algorithms 
used in LAPIN is given by Varner, et al. [2]. 

The governing equations used in the boundary layer 
code were the parabolized Navier-Stokes Equations 
(PNS). The PNS equations were 2-dimensional, vis- 
cous. unsteady, and were solved for the primitive vari- 
ables. The y-momentum equation simply reduced to 


x= 


x 


( 6 ) 


dP 


=0 


(16) 


The other governing equations for the boundary layer 
were as follows: 


2 


continuity 


3(p u) x 3(p v) + d p_ p 
dx dy dt 


(17) 


x-momentum 
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energy 
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The variables were defined as follows: 


tions, they were solved at each time step by marching 
from station to station in the direction of the core flow 
(x-direction). Solutions at a given x-station were ob- 
tained by solving the boundary layer equations sequenti- 
ally. First, the momentum equation was solved for the 
velocity component in the x-direction (u). Second, 
continuity was solved for the velocity component in the 
y direction (v). Then, the energy equation was solved 
for the temperature. If the solution had not converged, 
the momentum, continuity, and energy equations were 
solved again for the velocity components and tempera- 
ture, using the flow properties from the previous itera- 
tion. Convergence at a given x-station was obtained 
when the x-velocity components (u) at all y-grid loca- 
tions at that x-station converged. After the values at the 
station converged, the algorithm marched downstream to 
the next x-station. This streamwise marching continued 
throughout the entire solution domain. The differencing 
of the governing equations followed that of Kwon. el al. 
[3]. A detailed description of the differencing and 
solution algorithm used in the boundary layer algorithm 
is given by Roach, et al. [IJ. 


P= n 

Pmax 

(20) 

u~ u 

^max 

(21) 

P 2 

Pmax u max 

(22) 

7* C P T 

U m , K 

(23) 

C=t u^ x 

(24) 

Pinax^Wx 

(25) 




R 0 is the radius of the wall for axisymmetric flow (k=I). 
For rectangular coordinates R 0 was ignored <k=0). It 
should be noted that the Re and time terms were not 
actually non-dimensional. Re had units of (length)* 1 , 
while the time term had units of length. Also, the 
length variables (x, y. and Rj remained dimensional. 

So, each term in each of the above boundary layer equa- 
tions had units of (length) [ . 

The boundary layer equations were solved at a given 
time step knowing the solution at the previous time step. 
Due to the parabolic nature of the boundary layer equa- 


The boundary layer code used a modified Baldwin- 
Lomax turbulence model. Modifications were made for 
regions of adverse pressure gradients and separation. A 
description of the modified Baldwin-Lomax model is 
given in Sakowski. et al. [4]. 

In addition, downstream effects were propagated 
upstream in the subsonic portions of the boundary layer 
by using the pressure sweep and a weighted for- 
ward/backward differencing of the pressure gradient of 
Davis and Barnett [5]. 

III. Interface of Core Flow and Boundary Layer 

As described in the previous section, the core flow 
and boundary layer models dealt with the solution vari- 
ables in different manners. The core flow model used 
conservation forms, while the boundary layer model 
used non-conservation forms of die governing equations. 
Also, the two models non-dimensionalize their variables 
in very different manners. Thus, one of the critical 
steps in interfacing the two schemes is just translating 
values from one scheme to the other. The boundary 
layer model had to be able to interpret the core flow 
values, and. likewise, the core flow model had be able 
to interpret the results of the boundary layer model In 
addition to this bookkeeping, several other consider- 
ations were important. First, the interaction between the 
inviscid core and the viscous boundary layer had to be 
modeled quickly and accurately. Second, the coordinate 
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system had to be set up to maintain the parabolic nature 
of the boundary layer equations in the x-y plane, so that 
the algorithm could march in the x-direction. Third, the 
boundary layer properties had to approach the core flow 
values smoothly as we moved away from the wall of the 
inlet. Finally* die interface could not compromise the 
numerical stability of the combined viscous/inviscid 
model. 


The transform used is 

y T =y*arf{x, t) (26) 

the transformed x and t are the same as the original x 
and t. 


x T =x 


(27) 


Iterating between the core flow and boundary layer 
codes at a given times step was found to be time con- 
suming and unnecessary for the purpose of this model 
since such iteration had little effect on the overall flow 
behavior. A quick method of modeling the core 
flow/boundary layer interaction was found to be suffi- 
cient in this case. The core flow/boundary layer interac- 
tion was set up with a time staggering between LAPIN 
and the boundary layer codes. First, at a given time 
step, the predicted boundary layer parameters from the 
previous time step (displacement thickness and friction 
coefficient) were used to calculate the core flow proper- 
ties. Then the core flow properties were used to calcu- 
late the boundary layer properties at the next time step. 

The coordinate system for the boundary layer algo- 
rithm had to be set up to maintain the parabolic nature 
of the boundary layer equations in the x-y plane. One 
of the boundary layer assumptions was the velocity in 
the direction perpendicular to the core flow is much 
smaller than the velocity in the core flow direction. In 
many boundary layer algorithms, the coordinate system 
for the boundary layer was set up to follow the wall, 
with x along the wall and y perpendicular to the wail. 
There would be no problem with this approach if the 
duct had a constant area or if the core flow was solved 
as two dimensional. But, inlets were converging and 
diverging and LAPIN modeled the core flow as one 
dimensional. In an actual inlet the inviscid core flow 
far from the wall would not be parallel to the wall. 

Thus, as the two dimensional boundary layer code 
moved away from the walk the velocity perpendicular to 
the wall may not necessarily have been small compared 
with the velocity parallel to the wall. In order to match 
the one dimensional core flow and maintain the parabol- 
ic nature of the equations, a coordinate transformation 
was set up, based on the cove flow stream area. The 
coordinate system followed the stream lines of the l-D 
inviscid core flow. This way the velocity component 
perpendicular to the core flow stream direction was 
small everywhere, and the boundary layer assumptions 
could still be applied. 


The function arf(x,t) can be taken to be one over the 
core flow stream area. The boundary layer equations 
were solved in the transformed (x r . y r . t f } space. The 
problem with this coordinate system is that x r was not 
really perpendicular to y r . This comes from the fact 
that the core flow solution is l-D (x-direction only), 
while the real inlet would have had flow in both the x 
and y directions. This coordinate system was well 
suited for the convective terms; however, certain terms 
should have been expressed in the real (x, y, t) space, 
such as the viscous terms in the x-momentum equation. 
For example, the shear stress should have acted act 
along lines of constant y, not lines of constant yk For 
these terms we needed to transform back to the real 
space to get their values. To do this, the chain rule was 
used. The partial derivatives in the original boundary 
layer equations could be replaced with derivatives in the 
transformed space. 

For a generic scalar parameter (J), 



=arf{x T , t T ) -Q- 
dy T 


(29) 



d$ y 7 darf{x T ,t T ) {?) q ) 

dy 7 arf (x 7 , t T ) dx T 

There was some approximation associated with this 
scheme, such as mixing terms for the (x. y, y) space and 
the (x 1 , y\ t 1 ) space. Again, this is a consequence of 
matching a 2-D boundary layer model with a l-D core 
inviscid model of a flow that was not completely l-D. 
This transformation should not be necessary if we were 
using a 2-D core flow solution. Even with approxima- 
tion, this transformation worked well as long as there 
was no y-momentum equation. The v that was calculat- 
ed in the transformed space fit into the advective terms 
of the x-momentum and energy equations. However, a 
yLmomentum equation would be ill-posed since v 1 is 
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not always in the same direction nor perpendicular to u 
or u T . 


layer algorithm will approach the core flow values as 
the y-derivatives approach zero. 


Another concern when interfacing the boundary layer 
algorithm with the core flow algorithm was that the 
boundary layer properties had to smoothly approach the 
core flow values. That is, when the derivatives perpen- 
dicular to the wall were zero, the core flow values had 
to be a solution to the boundary layer equations. The 
complicated part of matching the boundary layer and 
core flows is the two ways the boundary layer and core 
flow algorithms were differenced. If the differential 
equations were solved exactly there would not have 
been a problem, but they were not solved exactly. LAP- 
IN used non-dimensional forms of p. pu. and p(e + l 2 u 2 ) 
as the solution variables, while the boundary layer algo- 
rithm used non-dimensional forms of the primitive vari- 
ables p. u, and T. Thus, what the core flow algorithm 
predicted as a solution, may not have been exactly what 
the boundary layer algorithm predicted as a solution as 
the y-derivatives go to zero (far from the wall). The 
difference is usually fairly small (2% or so), but this 
small difference can have a big effect on the integral 
performed to calculate the displacement thickness. 
Adjustment of how the pressure gradient term was cal- 
culated in the boundary layer algorithm forced the boun- 
dary layer properties to smoothly approach the core 
values. Without this adjustment, the displacement thick- 
ness, predicted by the boundary layer program had large 
errors. The adjustment of the pressure gradient was 
performed by solving the finite differenced momentum 
equations for dP/dx when all y-derivatives were zero. 
Similarly, the finite-differenced energy equation is 
solved for dP/dt when all y-derivative were zero. Also, 
without this adjustment, small errors in the velocity 
profiles at each station were additive as the algorithm 
marched downstream and marched in time, resulting in 
large errors in the friction coefficient near the exit plane 
of the solution domain and at later time steps. 


P„ dUifuJ 

~5x * ofE R k dx (31) 


JP-dT-^ p, u,d(R 0 *U- a ) _ap 
5t V * 51 R k dx 


(32) 


These corrected values of the pressure derivatives 
were calculated with the edge values from the core flow 
algorithm using the same differencing scheme used in 
the boundary layer algorithm. In this way the boundary 


The last consideration for interfacing LAPIN and the 
boundary layer code is a stability consideration. At 
times the boundary layer algorithm has a stability prob- 
lem. This problem tended to initiate near the edge of 
the boundary layer. From one iteration to the next the 
values near the core flow sometimes fluctuated between 
less than and greater than the core flow value. Some- 
times these fluctuations died out and the program con- 
verged. However, other times the oscillations grew, 
causing the calculations to diverge. To solve this prob- 
lem flaring was used. RVB was part of the advective 
terms in the y-direction. RVB was the coefficient in 
front of the <3u/<?y term in the finite differenced x-mo- 
mentum equation and the dT/dy term in the finite diff- 
erenced energy equation. These were the advective 
terms in the y-direction. Without flaring RVB was sim- 
ply pv. With flaring RVB was changed as follows: 


RVB=K X \ p v | 


(33) 


K,'s for momentum equation 

if u<u 9 if u>u 9 


if 





if 




*k — 


u 


K/s for energy equation 

if T<T e if t>t 9 
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The motivation for the above flaring was to make the 
core flow value a numerically stable solution in the 
boundary layer algorithm far from the wall. This flar- 
ing was found to be very important to help the stability 
of the algorithm, particularly when there was an adverse 
pressure gradient, separation, bleeds, or bypasses. 

IV. Results and Discussion 

Before describing the results of running the boundary 
layer/LAPIN codes for an actual inlet, we will look at 
two simple illustrative examples to show the effects of 
the corrections described in the last section. 

The first illustration is qualitative. Three simple cases 
were run. All three modeled steady flow over a flat 
plate. In the first case the correct pressure gradient was 
used. In the second case the pressure gradient was 
increased to 0.02/ft over the correct value. Finally, the 
third case was rim by decreasing the pressure gradient 
by 0.02/ft. Figure I shows the predicted velocity pro- 
files 2.5 feet upstream of the leading edge. At this 
point, the edge Mach number was 1.3. As shown on 
Figure 1. when the correct pressure gradient was used, 
the velocity profile smoothly approached the core flow 
value (u/u^=l). However, if the pressure gradient was 
too big, the velocity profile approached a value less than 
the core flow value as the y-derivatives went to zero, far 
from the wall (then jumped to the core flow value at the 
last point because of the specified boundary condition). 
Meanwhile, if the pressure gradient used was too small, 
the velocity profile approached values larger than the 
core flow values as the y-derivatives approached zero, 
far from the wall (then again jumped to the core flow 
value at the last point because of the forced boundary 
condition). 

This error in the pressure gradient impacted the pre- 
dicted friction coefficient. The friction coefficient pre- 
dicted when the pressure gradient was too low was 26% 
higher than the friction coefficient predicted when the 
pressure gradient was too high. However, the more 
significant impact was on the displacement thickness. 
When the pressure gradient was too high, the displace- 
ment thickness was 2.4 times the displacement thickness 
predicted when the correct pressure gradient was used. 
When the pressure gradient was too low. an unrealistic 
negative displacement thickness was predicted. The 
error in the displacement thickness was a strong function 
of the height (y-direction) of the solution domain. 

As a second illustration, a case was run for unsteady 
flow over a flat plate. In this illustration the boundary 


layer is not interactive with the core flow. Figure 2 
shows Mach number distribution at a particular time 
step. Two cases were run, one using the pressure gradi- 
ent from the core flow solution and one using the cor- 
rected pressure gradient calculated from equation 3 1 . 
Figure 3 shows the pressure gradients predicted by the 
core flow calculations and the corrected pressure gradi- 
ent from equation 32. The correction in the pressure 
gradient was small and had only a slight impact on the 
predicted Diction coefficient (figure 4). However, as 
shown on figure 5 the uncorrected pressure gradient 
allowed the boundary layer code to predict unrealistic 
negative displacement thicknesses when the pressure 
gradient was favorable. Then the displacement thick- 
ness predicted using the core flow pressure gradient 
diverged from the displacement thickness predicted 
using the corrected pressure gradient when the pressure 
gradient was adverse. 

In this second illustration, without the flaring de- 
scribed in the previous section the algorithm was un- 
stable. 

Now we will look at a more complicated example. 
LAPIN was combined interactively with the boundary/ 
layer code as described in the previous section. The 
combined model was used to predict flow in the NASA 
- Lewis 40-60 Inlet. The geometry of the inlet is given 
on Table 1 and pictured on figure 6, where r t is the 
radius of the cowl and t; is the radius pf the centerbody. 


Table I. NASA - LeRC 40-60 Inlet Geometry 


x/R t 

'/R, 

r/R c 

0.000 

0.000 


2.009 

0.446 

1.000 

2.152 

0.477 

1.000 

2.294 

0.509 

1.000 

2.437 

0.540 

1. 000 

2.580 

0.572 

0.997 

2.722 

0.604 

0.990 

2.865 

0.636 

0.983 

3.008 

0.665 

0.976 

3.150 

0.688 

0.970 

3.293 

0.705 

0.963 

3.436 

0.708 

0.954 

3.578 

0.698 

0.943 

3.721 

0.684 

0.933 
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3.864 

0.670 

0.923 

4.007 

0.656 

0.915 

4.150 

0.640 

0.908 

4.292 

0.623 

0.904 

4.434 

0.605 

0.902 

4.577 

0.586 

0.900 

4.720 

0.567 

0.900 

4.862 

0.546 

0.900 

5.005 

0.524 

0.902 

5.148 

0.501 

0.904 

5.291 

0.475 

0.906 

5.433 

0.449 

0.908 

5.576 

0.426 

0.912 

5.719 

0.403 

0.916 

5.861 

0.378 

0.918 

6.004 

0.351 

0.918 

6.147 

0.330 

0.918 

6.289 

0.309 

0.915 

6.432 

0.289 

0.908 

6.575 

0.271 

0.901 

6.717 

0.256 

0.893 

6.860 

0.245 

0.887 

7.003 

0.240 

0.882 

7.145 

0.239 

0.878 

7.288 

0.239 

0.873 

7.43 1 

0.239 

0.869 

7.573 

0.239 

0.865 

7.716 

0.239 

0.862 

Rc 

= 0.78 ft 


With bleeds at 

the following locations: 

3.223 < x/R c 

< 3.254 

on the cowl 

3.683 < x/R t 

< 3.697 

on the cowl 

3.783 < x/R c 

< 3.797 

on the cowl 

3.363 < x/R c 

< 3.429 

on the centerbody 

3.359 < x/R c 

< 3.373 

on the centerbody 

3.693 < x/R ( 

< 3.707 

on the centerbody 


And, a bypass at the following location: 

6.235 < x/R^ < 6.845 on the cowl 

The inlet was modeled with an upstream Mach num- 
ber (M*>) of 2.5 and with the exit Mach number chang- 
ing from 0.361 to 0.397 over 0.005 seconds. A time step 
of 0.0002 seconds was used. Figure 7 shows the pre- 
dicted displacement thickness after the tenth time step (t 
= 0.002 sec). The displacement thickness predicted 
without using the pressure correction is erratic and 
unrealistic. Because the boundary layer and core codes 
were interactive, the pressure correction now effects the 
combined algorithm as well. As shown on figure 8, the 
displacement thickness predicted without the pressure 
correction moved farther from the corrected model as 
the algorithm marched in time. In fact the algorithm 
actually became unstable on the 16th time step, so it did 
not predict any more solutions. 


As in the second illustrative example above, the algo- 
rithm was unstable without the flaring of the advective 
terms described in the previous section. 

Finally, the combined core flow/boundary layer model 
was compared with experimental results from Cubbison, 
et al. [6], run on the NASA - Lewis 40-60 inlet. The 
Mach number upstream of the inlet (ML) was again 2.5. 
Figure 9 shows the ratio of local static pressure to the 
free stream static pressure versus position in the inlet 
when a normal shock was near the throat of the inlet. 
Without the boundary layer code LAPIN predicted the 
shock wave to be downstream of experimental position. 
The combination of LAPIN and the boundary layer code 
predicted the shock wave farther upstream, closer to the 
experimental value. Also, the prediction of the exit 
pressure is improved when the boundary layer code was 
included. (Note: an exit Mach number boundary condi- 
tion was used). Figure 10 shows similar results when 
the shock wave was well downstream of the throat. 
Again, with the boundary layer code, the prediction of 
the shock position is closer to the experimental value 
than without the boundary layer code. Also, the predic- 
tion of the pressures near the throat was also improved 
by including the boundary layer code. 

V. Conclusion 

The uncoupled boundary-layer scheme was shown to 
be an effective method for correcting for viscous effects 
in a model of a supersonic inlet. This was effective 
when the overall flow variables are of concern such as 
in a controls type application. However, a more de- 
tailed model of shock-wave boundary layer interactions 
and other specific flow details should be used for de- 
tailed designs of inlets. 

Corrections to the model were made to maintain the 
parabolic nature of the boundary layer equations and to 
ensure the boundary layer properties smoothly approach- 
ed the core flow values. Also flaring was used to im- 
prove the stability of the algorithm. With these correc- 
tions the boundary layer algorithm corrected for bound- 
ary layer effects near the throat of an inlet and near a 
shock wave. 
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Figure 1 Predicted velocity profiles using pre- 

scribed pressure gradients. 



8 0 12 3 4 5 

X [ft] 


Figure 2 Core flow Mach number versus x 

for flow over a flat plate with a 
variable pressure gradient. 
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Figure 3 Corrected pressure gradient (equation 31) 
and uncorrected pressure gradient (using 
core flow values) for flow over a flat 
plate. 
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Figure 4 Predicted friction coefficient with and 
without the pressure gradient correction 
of equation 31. 
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Figure 5 Predicted displacement thickness 

with and without the pressure gradi- 
ent correction of equation 31. 
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Figure 6 Geometry of the NASA - LeRC 40- 

60 supersonic in lei. 
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Figure 7 Predicted displacement thickness on the Figure 9 
cowl inside the 40-60 inlet with an with- 
out the pressure gradient correction of 
equation 31, at time = 0.002 seconds. 


Predictions of pressure distribution 
within the 40-60 supersonic inlet, 
with and without the boundary layer 
algorithm, compared with experi- 
mental results of Cubbison, et ai. 

[ 6 ], 







*/R c 



2 3 4 5 6 7 8 


x/R c 


Figure 8 Predicted displacement thickness on the Figure 10 
cowl at location x/R c =329. 


Predictions of pressure distribution 
within the 40-60 supersonic inlet, 
with and without the boundary layer 
algorithm, compared with experi- 
mental results of Cubbison. et al. 
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